#ifndef _SIUT_SIMD_STANDARD_BBOX_HPP_
#define _SIUT_SIMD_STANDARD_BBOX_HPP_

#ifdef DEBUG_BBOX
#define THROW_BBOX
#include <strstream>
#include <stdexcept>
#endif

#include "Vec3f.hpp"


#ifndef FLT_EPSILON
#define FLT_EPSILON 0.00001
#endif

namespace siut
{
  namespace simd
  {

    struct BBox3f
    {
      simd::Vec3f min_;
      simd::Vec3f max_;
      bool empty_;

      BBox3f()
	:empty_(true)
      {
	;
      }
#ifndef NVCC
      BBox3f(simd::Vec3f min, simd::Vec3f max)
	:min_(min),
	 max_(max),
	 empty_(false)
      {
#ifdef THROW_BBOX
	sanityCheck();
#endif
      }     
#endif
    };
    ////////////////////////////////////////////////////////////////////////////////////////////////////
    //forward declaration
    inline simd::Vec3f closestPoint(const BBox3f &b, const simd::Vec3f &p);
    inline bool conjunct(const BBox3f &b, const BBox3f &bbox);
    inline bool containedClosed(const BBox3f &b, const BBox3f &bbox);
    inline bool containedClosed(const BBox3f &b, const simd::Vec3f &p);
    inline bool containedOpen(const BBox3f &b, const simd::Vec3f &p);
    inline float content(const BBox3f &b);
    inline simd::Vec3f corner(const BBox3f &b, size_t n, Vec3f &corner);
    inline float delta(const BBox3f &b, size_t i);
    inline float diagonal(const BBox3f &b);
    inline bool disjunct(const BBox3f &b, const BBox3f &bbox);
    inline float euclideanHausdorffDistance(const BBox3f &b, const BBox3f &bbox);
    inline simd::Vec3f getCenter(const BBox3f &b);
    inline const simd::Vec3f& getMax(const BBox3f &b);
    inline const simd::Vec3f& getMin(const BBox3f &b);
    inline simd::Vec3f getSizes(const BBox3f &b);
    inline void include(BBox3f &b, const BBox3f &bbox);
    inline void include(BBox3f &b, const simd::Vec3f &p);
    inline bool intersectBBoxes( const BBox3f &b, const BBox3f &bbox);
    inline BBox3f intersection(const BBox3f &b, const BBox3f &bbox);
    inline bool isEmpty(const BBox3f &b);
    inline simd::Vec3f lerp(const BBox3f &b, const simd::Vec3f &p);
    inline void sanityCheck(const BBox3f &b);
#ifdef DEBUG_BBOX
    std::ostream& operator<<(std::ostream &os, const BBox3f &b);
#endif
    ////////////////////////////////////////////////////////////////////////////////////////////////////

    inline simd::Vec3f closestPoint(const BBox3f &b, const simd::Vec3f &p)
    {
      Vec3f ret;
      for(size_t i = 0; i < 3; i++)
	{
	  if(p[i] < b.min_[i])
	    {
	      ret[i] = b.min_[i];
	    }
	  else if(p[i] > b.max_[i])
	    {
	      ret[i] = b.max_[i];
	    }
	  else
	    {
	      ret[i] = p[i];
	    }
	}
      return ret;
    }

    inline bool conjunct(const BBox3f &b, const BBox3f &bbox)
    {
      return !disjunct(b, bbox);
    }

    inline bool containedClosed(const BBox3f &b, const BBox3f &bbox) 
    {
      if(b.empty_)
	{
	  return false;
	}

      if(containedClosed(b, bbox.min_) && containedClosed(b, bbox.max_))
	{
	  return true;
	}
      else
	{
	  return false;
	}
    }

    inline bool containedClosed(const BBox3f &b, const simd::Vec3f &p)
    {
      if(b.empty_)
	{
	  return false;
	}

      for(size_t i = 0; i < 3; i++)
	{
	  if(p[i] + FLT_EPSILON < b.min_[i] || p[i] + FLT_EPSILON >= b.max_[i])
	    {
	      return false;
	    }
	}
      return true;
    }
     
    inline bool containedOpen(const BBox3f &b, const simd::Vec3f &p)
    {
      if(b.empty_)
	{
	  return false;
	}
      for(size_t i = 0; i < 3; i++)
	{
	  if(p[i] < b.min_[i] || p[i] >= b.max_[i])
	    {
	      return false;
	    }
	}
      return true;
    }

    inline float content(const BBox3f &b)
    {
      float ret = 1.f;
      for(size_t i = 0; i < 3; i++)
	{
	  ret *= (b.max_[i] - b.min_[i]);
	} 
      return ret;
    }

    inline simd::Vec3f corner(const BBox3f &b, size_t n, Vec3f &ret)
    {
      for(size_t i = 0; i < 3; i++)
	{
	  ret[i] = (n >> i)&1 ? b.min_[i] : b.max_[i];
	}
      return ret;
    }

    inline simd::Vec3f corner(const BBox3f &b, size_t n)
    {
      simd::Vec3f ret;
      for(size_t i = 0; i < 3; i++)
	{
	  ret[i] = (n >> i)&1 ? b.min_[i] : b.max_[i];
	}
      return ret;
    }

    inline float delta(const BBox3f &b, size_t i) 
    {
      return b.max_[i] - b.min_[i];
    }

    inline float diagonal(const BBox3f &b) 
    {
      return length(getSizes(b));
    }

    inline bool disjunct(const BBox3f &b, const BBox3f &bbox)
    {
      if(b.empty_ || bbox.empty_)
	{
	  return true;
	}
      for(size_t i = 0; i < 3; i++)
	{
	  if(b.min_[i] > bbox.max_[i] || bbox.min_[i] > b.max_[i])
	    {
	      return true;
	    }
	}
      return false;
    }

    inline float euclideanHausdorffDistance(const BBox3f &b, const BBox3f &bbox)
    {
      float dd, d = 0.0;
      simd::Vec3f t1, t2;
      size_t i;
      for(i = 0; i < 3; i++)
	{
	  t1 = corner(b, i);
	  t2 = closestPoint(b, t1);
	  dd = euclideanDistanceSquared(t1, t2);
	  if(dd > d)
	    {
	      d = dd;
	    }
	}
      for(i = 0; i < 3; i++)
	{
	  t1 = corner(bbox, i);
	  t2 = closestPoint(bbox, t1);
	  dd = euclideanDistanceSquared(t2, t1);
	  if(dd > d)
	    {
	      d = dd;
	    }
	}
      return std::sqrt(d);
    }

    inline simd::Vec3f getCenter(const BBox3f &b) 
    {
      return lerp(b.min_, b.max_, 0.5);
    }

    inline const simd::Vec3f& getMax(const BBox3f &b)
    {
      return b.max_;
    }

    inline const simd::Vec3f& getMin(const BBox3f &b)
    {
      return b.min_;
    }
     
    inline simd::Vec3f getSizes(const BBox3f &b) 
    {
      return b.max_ - b.min_;
    }

    inline void include(BBox3f &b, const BBox3f &bbox)
    {
      if(!bbox.empty_)
	{
	  include(b, bbox.min_);
	  include(b, bbox.max_);
	}
      else
	{
	  b.empty_ = false;
	  b.min_ = bbox.min_;
	  b.max_ = bbox.max_;
	}
    }

    inline void include(BBox3f &b, const simd::Vec3f &p)
    {
      if(b.empty_)
	{
	  b.min_ = b.max_ = p;
	  b.empty_ = false;
	}
      else
	{
	  for(size_t i = 0; i < 3; i++)
	    {
	      b.min_[i] = std::min(b.min_[i], p[i]);
	      b.max_[i] = std::max(b.max_[i], p[i]);
	    }
	}
    }

    /**
     * use conjunct instead
     */
    inline bool intersectBBoxes( const BBox3f &b, const BBox3f &bbox)
    {
      return !disjunct(b, bbox);
      /*
	if(b.empty_ || bbox.empty_)
	{
	return false;
	}
	for(size_t i = 0; i < 3; i++)
	{
	if(b.min_[i] > bbox.min_[i] || bbox.min_[i] > b.max_[i])
	{
	return false;
	}
	}
	return true;
      */
    }

    inline BBox3f intersection(const BBox3f &b, const BBox3f &bbox)
    {
      if(b.empty_ || bbox.empty_)
	{
	  return BBox3f();
	}

      simd::Vec3f bmin, bmax;
      bool nonempty = true;

      for(size_t i = 0; i < 3; i++)
	{
	  bmin[i] = std::max(b.min_[i], bbox.min_[i]);
	  bmax[i] = std::max(b.max_[i], bbox.max_[i]);
	  if(bmax[i] <= bmin[i])
	    {
	      nonempty = false;
	    }
	}
      if(nonempty)
	{
	  return BBox3f(bmin, bmax);
	}
      else
	{
	  return BBox3f();
	}
    }

    inline bool isEmpty(const BBox3f &b)
    {
      return b.empty_;
    }     

    inline simd::Vec3f lerp(const BBox3f &b, const simd::Vec3f &p)
    {
      simd::Vec3f ret;
      for(size_t i = 0; i < 3; i++)
	{
	  ret[i] = (1.f - p[i])*b.min_[i] + p[i]*b.max_[i];
	}
      return ret;
    }

#ifdef THROW_BBOX
    inline void sanityCheck(const BBox3f &b)
    {
      if(b.empty_)
	{
	  return;
	}
      for(size_t i = 0; i < 3; i++)
	{
	  if(b.min_[i] > b.max_[i])
	    {
	      stringstream s;
	      s << "min > max in " << __FILE__ << " at " << __LINE << std::endl;
	      throw runtime_error(s.str());
	    }
	}
    }
#endif

#ifdef DEBUG_BBOX
    inline std::ostream& operator<<(std::ostream &os, const BBox3f &b)
    {
      Vec3f corner;
      for (size_t i = 0; i < 3; i++)
	{
	  corner(b, i, corner);
	  os << corner << ' ';
	}
      return os;
    }
#endif

  }//end namespaces
}

#endif
